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SUMMARY 


An alternating direction implicit algorithm is presented for solving 
the conservative, full potential equation for unsteady, transonic flow. 

A new development is the time linearization of the density function. This 
linearization reduces the solution process from solving a system of two 
equations at each mesh point to solving just a single equation. Two 
sample cases are computed. First, a one dimensional traveling shock wave 
is computed and compared with the analytic solution. Second, a two 
dimensional case is computed of a flow field which results from a thicken- 
ing and subsequently thinning airfoil. The resulting flow field, which 
includes a traveling shock wave, is compared to the flow field obtained 
from the low frequency, small disturbance, transonic equation. 




TABLE OF SYMBOLS 


a 

speed of sound 

c 

characteristic length, e.g. chord length 

d 

characteristic time 

k = coc/q^o 

reduced frequency 
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pressure 

q> q 

pressure coefficient 
velocity, speed 
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Y 

specific heat ratio 
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velocity potential 
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gradient operator 

jump or difference in the quantity Q on the two sides 
of a surface of discontinuity 

Subscript 


00 

indicator of free stream values, e.g. M denotes free 
stream Mach number 
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I . INTRODUCTION 


Unsteady transonic flows are one of the most important yet least 
understood areas in fluid mechanics.^ For example, dips in flutter 
boundaries are often observed in the transonic Mach number range. An 
especially pronounced dip was observed for a supercritical wing.^ However, 
the mechanism for this dip is not yet understood. Other important problems 
in unsteady transonic flows are the flow about helicopter rotor blades,^ 
gust alleviation at transonic speeds, active control for dynamic stability 
at transonic speeds and transonic flow in turbines . 

This work is intended to provide a computational tool for the study 
of unsteady transonic flow and especially for the prediction of unsteady 
loads on aerodynamic bodies. We present an algorithm which solves the two 
dimensional conservative, full potential equation for unsteady flow. This 
equation will properly account for nonlinear effects, such as moving shock 
waves, under the assumptions that the viscous and rotational aspects of 
the flow are negligible. However it may be possible in the future to take 
into account those aspects in the case of a shock wave boundary layer 
interaction by coupling this algorithm with an unsteady boundary layer 
algorithm. 

This algorithm uses methods similar to those used to solve the low 
frequency, small disturbance equation.^”® As in that earlier development, 
certain criteria were used in choosing the methods. 

First, the algorithm should have fast computer run times, so that it 
can be routinely used for engineering applications. Hence the choice was 
made to develop a finite difference, alternating direction implicit 
algorithm. Such implicit algorithms permit large time steps; consequently 
the consideration of step size is based more on accuracy than on stability. 

Second, the algorithm should properly treat unsteady nonlinear effects, 
including moving shock waves. Presently, there is no well established 
method of shock fitting for unsteady flow. However, shock capturing 
methods are in common use.^’^ Therefore we choose the method of solving 
the governing equation in conservation form to ensure that the captured 
shock speed is the theoretically correct one. Calculations^ ’ based on 
equations in nonconservative form can produce captured shock speeds which 
depend on nonphysical parameters, such as mesh spacing or time step size. 

Third, the algorithm should properly handle aerodynamic shapes, 
arbitrary motions of those shapes, and arbitrary free stream Mach numbers. 

We choose the full potential equation as the governing equation. This 
equation removes three limitations^’”^ ’ ^ ^ on this third criterion which 
are assumed by low frequency, small disturbance, transonic, potential 
theory. These limitations are (1) that the flow field can be represented 
as a small perturbation about the free stream velocity, (2) that only low 
frequencies occur in the motion of the fluid and (3) that the free stream 
Mach number is close to one. This removal is especially necessary for the 
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accurate flutter analysis of supercritical airfoils. Such airfoils are 
thick and blunt nosed, and full potential theory properly handles these 
features . 

The fourth criterion considered was that the algorithm should be 
extendible to the three space dimensions. Alternating direction implicit 
methods are generally extendible 

Hence by the use of these criteria, we have decided to solve the 
unsteady, conservative, full potential equation by an alternating direction 
implicit method. This equation has two variables, viz. density and potential. 
So for a well posed problem we must augment the full potential equation with 
the unsteady Bernoulli equation to form a system of two equations. 

A new algorithm development is reported in this paper which simplifies 
the solution process for the resulting system of two finite difference 
equations at each mesh point. This development is the linearization back- 
ward in time of the nonlinear density function in the full potential 
equation. The linearization uses the unsteady Bernoulli equation and 
reduces the solution process to one which requires only the solution of 
a single equation at each mesh point for the potential. 

In Section II the assumptions made in deriving the governing equations 
are stated and those equations are presented in a general curvilinear 
coordinate system. In Section III the algorithm is tested on a one di- 
mensional traveling shock wave, for which the analytic solution is known. 

In Section IV a two dimensional flow is computed. Here a thickening 
and subsequently thinning airfoil produces a flow field with a rich 
structure, including a shock wave that travels in a large excursion. The 
computed results are found to be in close agreement with computed results 
from low frequency, small disturbance, transonic theory, as expected for 
this special case. 
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II. FORMULATION OF THE GOVERNING EQUATIONS 


The following assumptions are made on the flow of a perfect gas and 
the resulting nondimensionalized governing equations are given. 

1. Mass is conserved. 


k 


3t 


pq = 0 . 


( 1 ) 


2. The flow is irrotational . 

q = . 


( 2 ) 


3. The flow is isentropic. 

1 Y 
P = — 2 P 




(3) 


4. The flow is in viscid. 

From assumptions 2, 3 and 4, the momentum equation can be integrated^^ 
to yield the unsteady Bernoulli equation. 


3t 2 ^ Y-1 2 „2 


(4) 


M. (Y-1) 


Here the constant of integration is chosen so that for steady flow, Eq. (4) 
reduces to the steady Bernoulli equation. The variables have been non- 
dimensionalized in cartesian coordinates (t, x, y) as follows: x by c, 

y by c, t by d, <j) by q„c, a by q„, p by poo and p by p„q^. If there is a 
characteristic frequency co in the problem, then d = l/o) and k is the 
reduced frequency ; 


k = 


(5) 


Otherwise, k = 1, d = c/q^o and time is measured in characteristic lengths 
travelled at free stream speed q^. Here c is a characteristic length, 
which we choose as the chord length of the airfoil. 

Eq. (1) is in conservation form and hence has weak solutions, which 
conserve mass across surfaces where the flow variables are discontinuous. 
The jump condition, which is built into Eq. (1), is given by 

[kp] [pu] 11+ [pv] If = 0 , (6) 
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where the surface of discontinuity is given by S(x,y,t) = 0, and [ ] 
denotes the difference of the enclosed quantity on opposite sides of S. 

Computations of solutions to the flow equations about aerodynamic 
bodies are made more efficient by the use of a coordinate system that is 
adapted to the shape of the body and its motion. Hence we wish to write 
the flow equations, Eqs. (1) and (4) in a general coordinate system. Let 
be related to cartesian coordinates (t,x,y) by the transformation 

? = ? (x,y,t), n = n (x,y,t), t = t . (7) 

Then the conservation of mass equation in covariant form is given by 


k 




0 . 


( 8 ) 


See Appendix A for the derivation of Eq. (8). |j| is the absolute value 

of the Jacobian of the transformation (7), where 


J 






(9) 


U and V are the contravariant components of the velocity q in the ^ 
and q directions respectively, where 


U = 


+ 5 u + 5 V , 
X y 


V = rij. + n^u + HyV , 


( 10 ) 


and ^ = (u,v) in Cartesian coordinates. From Eq. (2), q = ('(•x» > 

hence by the chain rule for differentiation. 


u - > 


V - f , 


where , 


A, = 


2 2 

C + C , 

X y 


A., = 


C n + C h 
XX y y » 


A., = 


2.2 
ri + n 
X y 


( 11 ) 


( 12 ) 


5 


(13) 


Using Eqs. (3) and (4) and q 
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one obtains 


1 


P 





Equations (8) and (13) form the system of equations which will be solved 
in the domain of the variables (T,5,n). 

To complete the formulation of the problem, we specify the boundary 
and initial conditions. For inviscid flow, the boundary condition on an 
aerodynamic body is given by the physical requirement that the motion of 
the fluid is tangent to the body as the body moves. We will describe the 
body by 


of the velocity along n = 0 is zero; this condition is equivalent to 


In the sample computation of a two-dimensional nonlifting flow presented 
in this report, the far field boundary condition is that the flow variables 
have their steady free stream values. The initial conditions are obtained 
from the steady flow field, which is assumed to exist for tlO. 

Note that in a general flow problem, one can use an n equals con- 
stant surface to describe the location of an outer boundary. An example 
is flow in a wind tunnel, where an n -equals-constant surface locates a 
wall. Another advantage of using a general coordinate system is the 
ability to cluster mesh points in regions of interest, such as shock 
waves or stagnation points, without disrupting the smoothness of the mesh. 

The two dimensional calculations will be studied by plotting the 
resulting pressure coefficients, Cp , versus position on the airfoil. The 
formula for Cp is 


n (x,y,t) = 0 , 


(14) 


Mathematically, this boundary condition is that the physical component^® 


V = 0 . 


(15) 


Cp = (p^ - 1) . 


ym; 


(16) 


00 


Recall, before nondimens ionalizat ion Cp = (p - poo)/Q, where 
Q = (1/2) is the free stream dynamic pressure. 
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III. ONE DIMENSIONAL SHOCK WAVE MOTION 


First, we will compute the one dimensional case of a travelling shock 
wave. The shock wave travels downstream and separates supersonic flow 
upstream from subsonic flow. Here the state of the fluid is measured 
relative to the coordinate system and not relative to the moving shock. 

The exact analytic solution is known, so by comparison with the computed 
results, the shock capturing properties of the algorithm can be examined. 
In particular, we can verify that the algorithm is computing the correct 
shock speed, which is given by the jump condition, Eq. (6). 


a) Model Problem 

For this problem, we use cartesian coordinates, (x,t). Equations (8) 
and (13) simplify to 

+ (oKK = 0 ’ 

L XX 

The flow is along a finite interval and the length, L, of the interval is 
used for the characteristic length. Also, k = 1. The initial and boundary 
conditions are 

(0,t) = 0, (0,t) = U = 1, A (l,t) 

X Li X 

( X , 0 - X - xs , 

(|) (x,0) =j ° 

I xs + (x - xs ) Un , xs - X - 1, 
o o R o ’ 

xs^ = 0.045 and is the initial position of the shock. 

Here the velocity is nondimensionalized by the upstream (i.e. , left side) 
velocity; Figure 1 shows the initial conditions. 

The solution to Eqs. (17) and (18) is 

<|) (0,t) = 0 

0 < X < xs (t) , (19a) 

xs (t) < X < 1, 

xs (t) = xs^ + Vg • t, t > 0, (19b) 


(x.t) =< 


R 


Ur = 0.8. 


(18) 
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a 




Figure 1.— Initial Conditions for One Dimensional Problem. 


where xs(t) is the shock location at time t and Vs is the shock 
velocity. Vs is determined by the shock jump condition, Eq. (6), which 
for this case becomes 






V = 

S 




- 


( 20 ) 


where pl and are the upstream and downstream values of p. Using 
Eq. (17b) for pl and pR in Eq. (20), we see that Vs is determined 
by the parameters, u^, u^, ^‘*’t^R’ ^^t^R 


8 


are the upstream and downstream values of In our nondimensionalization 

of the variables in this one dimensional problem, we normalize by upstream 
values, so that u, = 1 and p = 1. From Eq. (17b), it follows that ((() ) = 0. 

Li Li t L 

In smooth regions of space time (x,t), the potential function (|)(x,t) 
satisfies 




( 21 ) 


'xt tx 

At a discontinuity, the shock jump condition which follows from Eq. (21) is 


fi 1 _ r . 1 3 S 

where S(x,t) = 0 is the curve of discontinuity in (x,t) . Since 

V ^ + 0 , 

s ax 3t 

u = 1 and (4> = 0, we have 

ti L Li 

<+t>R + \ <“r - 1) - 0- 

Using Eq. (24) in Eq. (17b), 

‘■e - + - 2 V 3(1 - u ^)}] '•-1 


( 22 ) 


(23) 


(24) 


(25) 


Using Eq. (25) in Eq. (20), we obtain a transcendental equation for Vg 
as a function of and u^. Since pl = Uj^ = 1, Eq. (20) simplifies to 


<»R - 1) ''s - 1 • (26) 

The following iterative scheme was used to determine Vg for the 
sample case computed, for which Moo = 1.2 and ur = 0.8. First, an initial 
guess of Tg was made. Next, Eq. (25) was used to compute pR. Then 
Eq. (26) was used to compute a new value of Vg. By a few choices of the 
initial guess for Vg , one was found which would generate a convergent 
sequence for Vs by iteration, using Eqs. (23) and (24) as described. 

The resultant Vq used was Vg = 0.4474296609271. This value of Vg is 
used in Eqs. (19) to complete the description of the analytic solution. 


b) 1-D Algorithm 

Next we present the algorithm used to compute this flow. First 
order accurate in time, implicit Euler, finite difference approximations 
to Eqs. (17) at the mesh point (i,n), located at the point (x,t) = 
(iAx,nAt) in the flow domain are 
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where 


n+1 

6^ p . +6 

t 1 ] 


/ n+1 


n+1 



(27a) 


(27b) 


(27c) 


6 = 


6 = 


6 = 
X 


6 = 
X 


(At)"^ (l - e“^) , 0 (At) 

(Ax)"^ (l - e"^) , 0 (Ax) 

(Ax)"^ (e^^ ■ ° 

(2Ax)"^ (e^^ ~ ^x^) ’ ° 


The shift operators are defined by 

+ 1 ,n 

*1 


r" ♦" 

X 1 


^n±l 

1 

n 


= h±i 


Eqs. (27) form a system of two implicit equations at each mesh point 
(i, n + 1) for the dependent variables arid A key idea 

developed in this report is the approximation of Eq. (27a) by an equation 
in which only the unknowns occur and they occur in a linear manner. 

The technique used for this approximation is the time linearization pro- 
cedure commonly used for the analysis and numerical solution of nonlinear 
differential equations. Although it is common practice^ to time 
linearize the nonlinear functions occurring in the space differences, 
here we will time linearize the nonlinear function occurring in the time 
difference in Eq. (27a) as well, by use of the special algebraic structure 
that the density possesses in Eq. (27b) (see acknowledgement). The resulting 
equation will preserve the property of conservation form, which is possessed 
by Eq. (27a) . 

First we will time linearize the nonlinear function in the space 
difference in Eq. (27a) by the following substitution for the density. 
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» Jti/2 ■ 01+1/2 * “ (“> • (28) 

Next we will time linearize the nonlinear function in the time difference 
in Eq. (27a) by the following Taylor series expansion for We have 


P i " Pi + 0 (At ) . 


Using Eq. (27b) for p^, we calculate (9p/3t)^ and substitute into the 
previous equation to obtain 


n+1 


n . 1 




(29) 


Substitution of Eqs. (28) and (29) into Eq. (27a) yields the first order 
accurate in time equation 







(30) 



(Pi+1/2 



= 0 . 


Note Eq. (30) is in conservation form, cf)j[ occur at level n+1 in a 
linear manner and do not occur at level n + 1. 


Using Eq. (27b) to determine p^, we now proceed to solve the scalar 
Eq. (30) at each mesh point. By algebraic operations, Eq. (30) is con- 
verted to the form 


(t)^'^^ + A6 + B6 (p^.t/,6 = C 

1 x^i x\^i+l/2 X / 


where 


~ n 

A = A t 6 (j) , 

X 


B . -4tV ' 


(31) 


K ) 

C - 2*'' - ♦”'2 + , 

D - {k ♦?)}] • 

Next we discuss our method of adding artificial dissipation. Numerical 
experiments indicate that Eq. (32) is a. stable finite difference scheme 
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for subsonic regions of flow. In supersonic regions we shall add an 
artificial dissipation term in order to damp numerical instabilities.^® 
The manner in which we add dissipation is to retard the evaluation of 
the density^ ^ ^ in the upwind direction, for the density term which 

occurs in the space difference. Hence in Eq. (31), Pi4-l/2 replaced 

by upwinded evaluation given by 


^ i+l/2 





i+l/2 


+ V . p 
1 


n 

i-l/2 • 


(32) 


The upwind weighting factor v? strongly affects the stability and 
accuracy of the computation. Insufficiently large values of can 

result in overshoots at shock waves or even unstable calculations. Too 
large values of result in solutions in which the shock profile is 

greatly smeared. Experience has shown that for unsteady flows, the 
proper values for depend on the computational mesh, the time step, 

the shock speed, and the shock strength. Currently the procedure for 
choosing is trial and error. In this example. 


= min jl, max 1 - 1 j c || , (33) 

where c is a positive constant, which for the case shown was tuned to 
7.5 and is the local Mach number. 



n\2 




Using Eq. (4), one obtains 



00 


(34) 


(35) 


In Eq. (33) the purpose of the min function is to retain an interpolation 
form in Eq. (32); otherwise unstable calculations can occur. For example, 
in the case presented here, removal of the min restriction caused the 
calculation to become unstable. 


Also, guided by the differencing of the low frequency, small dis- 
turbance, transonic equation®""^, certain terms were replaced by 
6x4) i terms in order to enhance the stability properties of the method 
The final algorithm is given by 


.n+1 , ^ .n+1 . ^ f 

+A6 (^. +B6 4). I 

1 XI x\ 1+1/2 X 1 / 


C 9 


(36) 


where 
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A = At S' <J>? , 

X 1 


n 


c = 2 (j) ; - (() . 


,n-l 





n 



+ M 





)^{ 


t d> 

X 1 


n-1 



c) Calculated Results 

The results of a calculation of the model problem are compared with 
the analytic solution in Figure 2. For the calculations, 101 equally 
spaced mesh points were used on the unit interval in the x direction. 
Also, since the velocity of the shock is known, the time step. At, was 
chosen so that the shock moves one-half mesh space per time step. The 
shock was initially located half-way between mesh points i = 5 and 
i = 6; xs(0) = 0.045. Equation (36) involves four levels of (j)±y so that 
three initial levels of data are required to start off the calculation. 
Equation (19b) was used to locate the shock at times xs (“At) and xs(“2At). 

The results are plotted for 49, 50, 99, and 100 iterations. Here 

(K)Li/ 2 = K-hl ' ’ 

and Ax = 0.01. The solid vertical lines show the shock position of the 
exact analytic answer. The results give numerical evidence that the 
algorithm is stable for Vg At 1 Ax/2, computes the correct shock speed 
and by the addition of the proper amount of artificial dissipation will 
produce a sharp shock profile. The ability to take a time step in which 
the shock moves one-half mesh spacing shows that the method allows the 
use of large time steps, by which we mean on the order of time steps 
allowed by the implicit algorithm of Ref. 6 for the low frequency, small 
disturbance equation. No determination of a limit on the time step size 
was made. 
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• n = 49 
□ n = 50 



Figure 2.— Shock Profiles for One Dimensional Problem. 
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IV. TWO DIMENSIONAL UNSTEADY TRANSONIC FLOW 


In this section we will present the results of a sample two dimen- 
sional calculation. The case presented will have no lift. However a 
strong shock wave will develop and subsequently travel rapidly upstream 
in a large excursion. Meanwhile the flow will become completely sub- 
sonic. The results will be compared with computed results from low 
frequency, small disturbance, transonic theory.^ 


a) Problem Formulation 

We will compute the flow that results from a thickening and sub- 
sequently thinning airfoil. Here time is measured in chord lengths 
travelled, k = l, and x and y are nondimensionalized by the chord 
length. A parabolic arc airfoil thickens from zero to 10 percent 
thick as a fluid particle travels 15 chord lengths, relative to the 
airfoil, at the free stream velocity q^. The airfoil then thins to 
zero thickness at 30 chord lengths of travel. Figure 3 shows the mid- 
chord thickness function 6(t), given by 


0.1 1 10- 15 (i) +6(5^) 'I , 0 i t < 15 , 


S(t) = 


0.1 
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+ 6(^)1 (^f . 15 < t < 30 , 


t > 30 . 


For our problem we can choose a simple, sheared coordinate system 
to fit the airfoil and its motion. Equations (7) become 

C = X , ri = y - S(x,t) , T = t , 


where 


y = S(x,t) , 0 ^ X ^ 1 , 


(37) 


(38) 


is the equation of the airfoil upper surface boundary in (x, y, t) 
space. S(x,t) 5 0 for x < 0 and x > 1. Since the flow is symmetrical 
about y = 0, we restrict the domain of flow to q > 0. The governing 
equations simplify. J = 1, so the conservation of mass equation, 

Eq. (1), becomes 


+ ^ “ 


(39) 


15 



Figure 3. — Midchord Thickness Function 6(t). 


where, using Eqs. (10) and (11) 


U = u 

V = -S - S u + V 

t X 


6 p — S 6 
^5 x^n 




■^t ■ = -S^ - + (1 + 


(40) 


The last equalities in Eqs. (40) follow from Eqs. (38); S(x,t) = S(? , x) , 
so and S^. = S^. The unsteady Bermoulli equation, Eq. (4), 

becomes 


P = 


1 + 
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The initial condition is q = for t 5 0, The body boundary con- 
dition, Eq. (15), becomes 


or 


V = + S u , 

t X ^ 


(1 + S^) 


( 42 ) 


Off the airfoil for n = 0, cj) =0 from symmetry. In the far field 
q = q^ = (u,v) = (1,0). From^Eqs. (40) 


u=(()j.— Scj) =1. 


(43) 


From Eqs. (38), (|) = c|)y, so (|) = c(>y = v = 0. Substituting into 

Eq. (43), we have^ (j)^ = 1 in ?he far field. 

For our computations, we approximate the flow on the infinite 
domain with the flow on a finite domain, 

-30 < 5 ^ 30 , 0 < n ^ 50 . (44). 


On the far boundaries, n = 50 and g = ± 30, we impose Dirichlet con- 
ditions, viz. (\> is specified with the result that = 0 on 
5 = ± 30 and ((>w = 1 on n = 50. The problem formulation is now 
complete. 


b) Implicit Factorization Algorithm 


The 2-D algorithm is a generalization of the ID differencing. 

Let (i,j,n) denote the indices of a mesh point located at the coordi- 
nates (5,n,x) in the flow domain. Here we use a stretched computa- 
tional mesh; see Fig. 4 for a sketch of the mesh. First order time 
accurate, implicit Euler, finite difference approximations of Eqs. (39) 
to (41) are 


n+1 , 


( ^ 


n+l 


^ n+l „n+l . A^+l 

*nW2,j 


+ 6 

n 


n+l 

^i,j+l/2 



6 d) 
n 


n+l 

i.j 



<P 


n+l 

l,j+l/2 



(45a) 
= 0 , 
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Figure 4.— Sketch of Sheared Computational Mesh. 



(45b) 


n+1 1/ n+1 

Pi+l/2,j " 2ri+l,j 


n+1 \ 
+ P . . 
1.1 / 


n+1 - i/ 

Pi,j+l/2 2\Pi,j+l 


+ P 


n+1 


). 


(45c) 
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where 


1 / .n+1 ^ .n+1 \ ^n+1 


(45d) 


= (^)r . 

(45e) 

= (At)"^ (1 - E"^) 

0(At) 



0(A5) 


^ «1+1 - “ 

0(AC) 


' <■ 4 ’’ - "f) • 

0(A^2) 



6 = 
T 

*5 


\ (1 - 


®n ' <"j+X - - « 


^ = (n . . 1 - n . , ) ^ (e"^^ - e ^) 

n j+1 j-1 n n 


0(An) 

0(An) 

0(An2) 


The shift operators are defined by 


.±1 


. n±l 

T 



,±1 




4> • • 

~ ^i±l,j 

,±1 


j.n 

'n 




Equations (45a) and (45b) form a system of two implicit equations 
at each mesh point (i,i, n+1) for the unknowns and . As in 

the one dimensional problem we linearize backward in time the density 
functions in Eq. (45a), so that only the unknowns occur and they 

occur in a linear manner. Again by a Taylor series expansion we have 


for that 

I5J 


Pm ■ pj.j * ml , “ ^ • 
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and substitute into 


Using Eq. (45b) for , we calculate (9p/3t)? , 

the previous equation to obtain 


C: = <j)- Cj) 

*ij)Nh <j)| 


1 „n+l 

2 


*li) 

- O 

<j) 

+ CK ♦lj)K O} +»<''''> • 


Also, in Eq. (45a) let 


n+1 n , n/ A \ 

^i+l/2,j ■ ^±+ 1 / 2,3 ^ 


n+1 


__ n 

'i,j+l/2 " ^i,j+l/2 


+ 0(At) 


(46) 


(47) 


As in the one^dimensional case, the 4>i,j terms in Eq. (46) are 
replaced by 6^ ^i, j terms, in order to enhance the stability of the 
difference scheme. 

The manner of adding artificial dissipation is similar to that in 
the one dimensional case. Again, we add an artificial dissipation 
term^® by retarding upwind the density evaluation. For this 
sample case, only densities occurring in the 5 difference term, viz, 
the second term in Eq. (45a), were retarded in the upwind direction. 
However, in general, one can retard densities upwind in both space 
difference terms in the manner shown in reference 19 and these achieve 
a rotated difference^® scheme. Here p?,i/o • is replaced by • > 

an upwinded evaluation given by 
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where 


'-n /- n \ n , n n 

Pi+l/2,j ~ (^ ■ '’i.jj Pi+l/2,j '"i.j ^±- 1 / 2,2 


( 48 ) 


V, . = min 
i.J 


1, max^ 


(49) 


where c is a positive constant, which for the computation shown was 
turned to 2.0 and . is the local Mach number. Recall M = q/a. 


so 


«.,) ' = ^ (\ *li) * ('n *U 


(50) 


(51) 


and using Eq. 4 


Rj)' ^ - < \*li) -Kof}- 


The (5,n,T) coordinate system is sheared and therefore nonorthogonal. 
Hence in Eq. (39), there are two cross derivative terms, viz ”(Sr (f)^)^ 
and “(Sg 4>^)q • In an implicit factorization algorithm, the question 
arises as to how to handle these terms. One would like to difference 
them implicitly on the assumption that such differencing would enhance 
the stability properties of the algorithm. However, here we will 
difference them explicitly. Recall that in alternating direction 
implicit schemes, each implicit difference equation is replaced by 
intermediate equations, which are typically implicit only along a 
single coordinate line. By differencing the cross derivative terms 
explicitly, the resulting intermediate equations in our algorithm will 
also link unknowns only on the same coordinate line. Hence the 

resulting algorithm will retain the structure of being completely 
vectorizable. This feature is desirable for the use of an algorithm 
on computers which are architecturally vector array processors. Hence 
in Eq. (45a) let 



Oi/al) ■ 


(53) 
+ 0(Ax) . 


Substituting Eq. 


(46) to (48) in Eq. (45a), we have 
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c ) - ♦ m )- *z] 

- K i\ ih oi - €'(\ *h } I'. (vm)I 

( 54 ) 


5^i,j 5 n^i+l/2,j 


) 


■" hi^l-H/2,j 

KKj^in |(^ + <^r">Vn*t3 - ^r«5*l3+l/2 - ®f'|) ■ 0 


Tl^”l 

Note Eq. (54) is in conservation form and (j)^ ^ occurs in a linear 

manner. We now shall factor Eq. (54). Note Eq. (45b) is used to 
evaluate the density at the previous time levels. 

By algebraic operations, Eq. (54) is transformed into the form: 


i,j ^^±,2 n "^1,3 5ri+l/2,3 ^^x,2l 


. n+1 


n+1 


. n+1 1 




(55) 


R , 


where 

A = A1 + A2 






B = B1 + B2 


B1 = Ax 


1 + S 


n+1 


6 ({>“ . 


„n+l t- , n 

*sh,j 


5 n+1 
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B2 = 6 (|)? i 6 (|)^ .] - , 

\ Cx \ C n ^ 1 , 2 ) 2 5 ) XT ’ 


C = -At 




2-Y 


A = A1 + A2 


A1 = E~^(A1), A2 = E"^(A2), B1 = E"^(B1) , 


B2 = E~^(B2), C = e;^^(C) , 


R = 2 (j)" . - + (Al) 6 ^ (|)? . + (Bl) 6 d)" . 

i,: i.J n ^x,j 

, 1 ) Ax / n n-l\ / n-l\ 2- y 

*“('") M's ♦l,j) 




2ip^ } 


+ A2 Sj + B2 6^ I j + C j sf" 6^ j) 


,n+l 


+ C 6 ( 

n v±. 


n+1/2 
3+1/2 


„n+l r ,n _n+l 

'c*i,J+l/2+^ 


)• 


Equation (55) is an implicit equation that determines (j) 


We wish to factor the operator L? . to obtain two intermediate 

i.J 

equations at each mesh point. Each intermediate set of equations 
will be implicit in only one coordinate direction. The resulting 
matrix equations will be tridiagonal and can be easily inverted by 
use of the Thomas algorithm. The factorization of the operators , 
will be approximate; however, the resulting additional terms will 
only be error terms that are of higher order error in time than the 
discretization errors already made. Also, some of the error terms will 
not be in conservation form. Therefore additional terms will be added at 
the nth time level to recover the conservation form for the factored 
equation. 


n+1 

i,j * 
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The factorization of Eq. (55) is 


5 implicit sweep 


= R - + (sf')')\ *",j) + CD , 


(56a) 


n implicit sweep 


' *n C <= (l ^ (^r)') % Ci ) 

■ *li + ®«n*i.3 * • 


where the term CD is added to the right side in the C sweep in 
order to maintain conservation form and 


(56b) 


D=^6. 


B6 5 (J)^ . + C6 , (l + (Sp)^)? 6 .) 

n T^i,j n\ i,j+l/2\ / n x"^i,j/ 


(56c) 


Note that the 5 implicit sweep is not a consistent approximation to 

Eq. (39) because of the replacement of the term B6^ c()^ . in Eq. (55) by 

n ^ 

the term B6^ (|)^ ^ in Eq. (56a). Equations (56) are a generalization of 

the alternating direction implicit (ADI) algorithm of Lees.^^ Equation 
(56b) can be solved for (j)J .. The resulting expression can be substi- 

j 

tuted into Eq. (56a) to yield 

h,3 ^ ^hY±+l/2,ihh.3l 


+ C 6 




(57) 




= R + CD . 
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Equation (57) is an approximation to Eq. (55); there are two ad- 
ditional terms on the left side, which resulted from the splitting and 
there is the additional term CD on the right side. Division Eq. (57) 
by the factor C and rearrangement of the terms results in the equation: 

T1 + T2 + T3 = 0 , (58) 

where T1 is all the terms on the left side of Eq. (54), 


T2 = At6 


n 


C Pi+l/2,j "5 


B 6 6 

n T "^i,j 




and 


(60) 


T3 - 6t- 6^ 


Using Eq. (56c) 


B 6 6 + C 6 (p? . , /, / 1 + 5 6 

q T ^1,3 n y*"!, 3+1/2 \ \C / / n t^i,3/ 


- D. 


(61) 


13 - '=«n(‘>i.3+l/2 (^rT)V'<5)l} . 


From the difference form of the terms Tl, T2 , T3, given in Eqs . (54), 

(59) and (61), we see that Eq. (58) is in conservation form, including the 
error terms, T2 and T3. Solutions of Eq. (57) are identical to solutions 
of Eq. (58), since only arithmetic operations are required to go from 
one equation to the other. Hence solutions of Eq. (57) will satisfy the 
conservative form, finite difference approximation, Eq. (58), to the con- 
servation of mass differential equation. 


Since A = 0 (At), B = 0 (At) and C = 0 (At^) , it follows from Eqs. 

(59) and (61) that T2 = 0 (At2) and T3 = 0 (At2) , So the additional 
error terms in Eq. (58) do not affect the first order accuracy in time of 
the finite difference approximation. 

The addition of the term CD in Eq. (56a) is required because the 
term At A6 c()^ . in that equation produces error terms which are in 
nonconservafive rorm. One alternative, which would eliminate that term 
from Eq. (56a) , is to replace the term 5g(()^^“l in Eq. (46) by j • 

However in the problem under consideration, tnat replacement resulted 
in an unstable algorithm. Another alternative is to interchange the order 
of the sweep directions in Eqs. (56) and also to replace 6^^ 

6^ (j>n . in Eq. (46). The rationale is that the 5 coordinate direction 
is app^roximately the flow direction in the supersonic region about an air- 
foil. For stability, in such regions, the j term must be differenced 
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implicitly. However the 6^ ^i,j term might be differenced explicitly 
and yet the algorithm would remain stable. The change in order of the C 
and n sweeps results in ?j>i/2,3 appearing in both sweeps. Hence the 
test for supersonic flow, Eq. (49), must be applied in both sweeps with 
resulting additional computations. However, in a rotated difference scheme,^® 
the density j+1/2 appearing in the difference term would also be 

upwind biased, so that interchanging the order of the sweeps would not 
introduce any additional testing. This second alternative has not yet 
been tested. 

c) Calculated Results 

The calculations were performed by solving Eqs. (45b) and (56) on 
the computational mesh. For the case shown, the mesh had 151 points in 
the 5 direction and 41 points in the n direction. There were 101 
mesh points in the 5 direction between 5 = -1 and 5 = +1 at a uniform 
spacing of = 0.02. The remainder of the K niesh on either end was 

smoothly exponentially stretched to the boundaries at 5 = ±30 chord 
lengths. The mesh in the n direction was smoothly, exponentially stretched 
to the far field boundary at n = 50 chord lengths. At j = 7, nj = 0.18 
and Apj = (hj-t-1 - h3-l)/2 = 0.05. So this mesh is very similar to the 
mesh used in Ref. 6 for the low frequency, small disturbance, transonic 
equation calculations. See Appendix B for a list of the computational 
mesh. 


As in the one dimensional case, three levels of initial data are 
required to advance the potential in time. The initial values i 

were determined by setting (|). ^ = 0 and requiring 6r ^f>i 1. Then we 

set ^i^j ” ^i,j ~ ^i j* boundary condition on the airfoil, Eq. (42) 

was imposed by the first order accurate approximation. 




n+l 

i,l 


T g K 


1 + 




n 

ill 


(62) 


where j = 1 mesh points lie on the coordinate line on the airfoil. Note 

that Eq. (62) is explicit in Off the airfoil, (4>p)5^]; = 0 was 

applied. At the far field boundaries, i = 1, i = 151 or j =41, 4>i,j was 
not changed from the initial values, with the result that j = 0 on 

i = 1 and i = 151 and . = 1 on j = 41 for all values of n*. In our 

calculations M =0.85 and’ A = 0.0625. 

00 X 

The pressure coefficients resulting from the calculation are plotted 
in Figs. 5 and 6. The airfoil is thickening for t < 15, as shown in Fig. 3, 
and in Fig. 5 we see that the effect of that thickening is the initial 

formation and downstream propagation of the shock wave. By t = 8.5, the 

flow has a small region of supersonic flow but there is no shock wave yet. 

By T = 11.5 a shock wave has formed; it is growing stronger and is moving 
downstream as the airfoil continues to thicken. At x = 18.25 the shock 
has slowed down and is nearly stationary. The local Mach number upstream 
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Full Potential Equation 
At = 0.0625 

Low Frequency, Small 
Disturbance, Transonic 
Equation, Ref. 6 
At = 0.125 

Moo “ 0.85 


Figure 5.— Plot of Pressure Coefficients for Formulation 
and Downstream Propagation of the Shock Wave. 
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Figure 6.— Plot of Pressure Coefficients for Upstream Propagation of the Shock Wave. 


of the shock, measured relative to the airfoil frame of reference, is 
approximately 1,38. Notice that the shock profiles are captured sharply 
by the algorithm. There was a time lag of about At = 3.25 between the 
time of the maximum thickness of the airfoil at t = 15 and the time of 
the shock strength reaching its maximum. 

The subsequent upstream movement of the shock wave is shown in 
Fig. 6. At T = 26.875 the flow has just become completely subsonic 
and the shock is now traveling rapidly upstream. The rapid excursion 
is shown by the plots at t = 29 and t = 32. In Fig. 6, the shock wave 
is from subsonic to subsonic flow and its strength is very weak compared 
to the shocks in Fig. 5. Since the computation of the density is no 
longer upwind biased during this later completely subsonic flow, the 
result is that the algorithm increasingly smears the shock profile. A 
calculation with a time step reduced to At = 0.03125 sharpened these 
weak shock profiles at t = 29 and 32. They appeared similar in shape 
to the low frequency, small disturbance profiles plotted in Fig. 6 but 
they were positioned near the centers of the full potential shock pro- 
files. At T = 11.5 and 18.25, the smaller time step calculation showed 
minor changes in the Cp profiles from the larger time step full potential 
calculation. The shock positions were moved slightly rearward and the 
shock strengths were slightly increased for the smaller time step case. 

Another view of the solution is given in Fig. 7, which is a plot of 
the midchord pressure coefficient as a function of time. There is a lag 
of about At = 2.75 between the maximum airfoil thickness and the 

maximum flow expansion time. The sharp drop in the curve occurs as the 
shock wave propagates upstream past the midchord point at about t = 26.6. 

There are no published results on unsteady, conservative, full potential, 
transonic calculations with which we can compare. Therefore we have com- 
pared our results with low frequency, small disturbance, transonic theory.^ 

The purpose of the comparison is twofold. First, for this problem, the 
two theories should be in close agreement. By comparison in Figs. 5, 6 
and 7, we conclude that they are in close agreement and therefore that 
there are no gross errors in the full potential calculations. 

The second purpose of the comparison is to see how the results of 
the two theories differ. Notice that at times t = 11.5 and t = 18.25, 
the low frequency, small disturbance, transonic shocks are rearward of 
the full potential shocks, whereas at t = 26.875 they are in close 
agreement. Solutions in small disturbance, transonic theory can be adjusted 
in various ways, such as by Krupp, Von Karman-Spreiter, or Cp* 
scaling. In our comparisons, Cp* scaling was used in the small 

disturbance calculations. In that scaling, the small disturbance critical 
pressure coefficient, Cp*, is equal to the steady full potential critical 
pressure coefficient. The result is that for steady flows that are barely 
supercritical, the resulting shock strengths and positions should be in 
close agreement. In Ref. 25, a plot shows a comparison of shock positions 
versus airfoil thickness for steady flows at a fixed Mach number between 
full potential theory and small disturbance theory. There Cp* scaling 


29 


U) 

o 


Full Potential Equation 
Ar * 0.0625 



Figure 7.— Plot of Midchord Pressure Coefficient Versus Time. 


was also used. That plot shows that the small disturbance shock position 
is rearward of the full potential shock position and that the separation 
increases as the thickness increases. Our comparison at x = 11.5 and 18.25 
shows the same relative position of the shocks. Also there is an increase 
in the separation of the two shocks at the later time, which can be 
attributed to a greater thickness effect if one allows for the time lag 
between the motion of the airfoil and the reaction of the fluid. At 
T = 26.875, the small disturbance flow is barely supercritical and the 
shocks have essentially the same strength and position, which is what 
would be expected from using Cp* scaling in a steady flow. 

A more detailed comparison of the two theories is made in Appendix C. 
There the effects on the pressure coefficient plots are shown for various 
approximations made to full potential theory in going to small disturbance 
theory. The approximations studied are modelling the airfoil as a slit, 
linearizing the boundary condition and linearizing the formula for Cp. 


31 


V. CONCLUDING REMARKS 


An algorithm has been developed which solves the two dimensional, 
conservative, full potential equation for unsteady, transonic flow by a 
finite-difference, alternating-direction implicit method. This develop- 
ment was motivated by the desire to remove three limitations of low- 
frequency, small-disturbance, transonic theory while retaining the 
computational efficiency and shock-capturing features of algorithms 
developed for that equation. This removal is necessary for an accurate 
treatment of thick, blunt nosed supercritical airfoils. A noteworthy new 
feature of the algorithm is the reduction of the system of two equations 
to be solved at each mesh point for the velocity potential and density 
functions to one equation involving only the unknown velocity potential. 

By first testing the algorithm on a one-dimensional traveling shock 
wave, we found, as shown in Figure 2, that the algorithm is stable, 
captures the weak solution shock motion correctly and allows relatively 
large time steps, similar in size to those permitted by the small- 
disturbance algorithm. Next we calculated a two-dimensional flow that 
was produced by a thickening and thinning airfoil. That test case, as 
shown in Figures 5-7, also is stable, allows large time steps, and is 
in close agreement with low-frequency, small-disturbance transonic 
equation calculations as required for this test case. 

Further work is underway. (1) A pilot code for use in the aero- 
dynamic community, is being prepared which treats arbitrarily shaped 
airfoils using an efficient body fitted mesh^^ and which allows for pitch- 
ing and plunging airfoil motions. This code will be capable of providing 
the unsteady loads for flutter and other aeroelastic calculations on both 
conventional and supercritical airfoils. (2) The algorithm is being 
extended to three space dimensions. (3) A modification to the algorithm 
is being investigated which would make the body boundary condition 
computation implicit. This modification could permit the use of even 
larger time steps. 
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APPENDIX A 


DERIVATION OF THE GOVERNING EQUATION IN A GENERAL 
CURVILINEAR COORDINATE SYSTEM 


The conservation of mass equation, Eq. (1) is written in cartesian 
coordinates ( t , x, y) as : 


9(kp) 3 (pu) 3 (pv) ^ ^ 

3t 8x 8y 


(Al) 


where q = (u,v) is the velocity vector in (x,y) space. We shall use 
tensor analysis^^'*^® to write this equation in a general coordinate 
system. In Euclidean space time, the line element ds^ in the Euclidean 
coordinates (t,x,y) is 


2 2 2 2 

ds = dt + dx + dy . (A2) 

To simplify the expressions that follow, we introduce tensor notation, 
for use in this appendix only. The coordinates of a point in space time 
will be denoted xV^, where the superscript y takes values 1, 2 or 3. 
Hence (t,x,y) = (x^,x^,x^). In tensor notation the line element ds^ 
has the form 


ds^ = g dx“dx^ , (A3) 

ap 

where here is the Kronecker delta function, and 

= 1 if a = B and = 0 if In Eq. (A3) we have used the 

Einstein summation convention, which assumes summation on "an index 
which appears twice in the general term."^® Hence in Eq . (A3) we sum on 
both a and B through the values 1 to 3. g<^g are the components of 
the metric tensor. 

Now consider the general transformation of coordinates 

T = x(t,x,y), 5 = 5(t,x,y), n = n(t,x,y) . (A4) 

Let x^ denote the coordinates of a point in this coordinate system, so 
(xl,x2,x3) = (T,?,n) and Eq. (A4) simplifies to x^ = x^(x"^). Then the 
line element ds^ in the coordinate system is 

ds^ = g dx^ dx^ , (A5) 

^yv 

Eq. (A5) follows from the use of the inverse transformation x^ = x^(x^^) 
so that 

djc" . ^ dx" . <A6) 

3x^ 
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Using Eq. (A6) in Eq. (A3) gives Eq. (A5) with 





9x' 


6 


r, a . a 
= ^ 9x 3x 

zL 3xVA 


9x^ 3x^ 


3xVA 9xv 


a=l 


(A7) 


Now we shall write the conservation of mass equation in a general coordi- 
nate system. Using the velocity vector q = (u,v) we form the contra- 
varient components of the velocity vector in space time in Euclidean 

coordinates (t,x,y) . 

= (k,u,v) . (A8) 

Also the gradient operator in (t,x,y) space time is defined by 

9 = (9,, 3,9). (A9) 

V t’ x’ y 

So in tensor notation, Eq. (Al) becomes 

9^ (PV*") = 0 . (AlO) 

Under the general transformation x^ = x^ (x^^) , Eq. (AlO) takes the co- 
varient form 

(^) . 0 . (All) 

where |j| is the absolute value of the Jacobian of the transformation. 


J = 



(A12) 


3x^/3x^ is the transformation matrix. are the contravarient components 

of the space time velocity vector in the xU coordinate system. 


3x^ 


(A13) 


Note that Eq. (All) is in conservation form. This form follows from an 
application of the fundamental Ricci theorem^® to the equation that results 
when the transformation of coordinates, Eq. (A4) , is applied to Eq. (AlO). 
The Ricci theorem states that the covarient derivative of the metric tensor 
is identically zero. 
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If T(t,x,y) =t, then the following simplifications result. 


j = C n - C n 
X y y X 


= k 


U = + CyV , 

-3 

v=v =n 4-nu + nv, 

L X y 


and Eq. (All) becomes 


(A14) 


(A15) 


Eqs . (A14) and (A15) are the desired equations and are used in Section II. 
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APPENDIX B 


COMPUTATIONAL MESH 


a) 


5 — Direction Mesh Distribution 


i I 

1 -.300000 E+02 

2 -.237062 E+02 

3 -.187520 E+02 

4 -.148554 E+02 

5 -.117937 E+02 

6 -.939063 E+01 

7 -.750718 E+01 

8 -.603333 E+01 

9 -.488218 E+01 

10 -.398506 E+01 

1 1 -.328770 E+01 

12 -.274720 E+01 

13 -.232969 E+01 

14 -.200837 E+01 

15 -.176209 E+01 

16 -.157411 E+01 

17 -.143122 E+01 

18 -.132297 E+01 

19 -.124112 E+01 

20 -.1 17913 E+01 

21 -.113186 E+01 

22 -.109522 E+01 

23 -.106601 E+01 

24 -.104167 E+01 

25 -.102020 E+01 

26 -.100000 E+01 

27 -.980000 E+00 

28 -.960000 E+00 

29 -.940000 E+00 

30 -.920000 E+00 

31 -.900000 E+00 

32 -.880000 E+00 

33 -.860000 E+00 

34 -.840000 E+00 

35 -.820000 E+00 

36 -.800000 E+00 

37 -.780000 E+00 

38 -.760000 E+00 

39 -.740000 E+00 

40 -.720000 E+00 

41 -.700000 E+00 

42 -.680000 E+00 

43 -.660000 E+00 

44 -.640000 E+00 

45 -.620000 E+00 

46 -.600000 E+00 

47 -.580000 E+00 

48 -.560000 E+00 

49 -.540000 E+00 

50 -.520000 E+00 


i ^ 

51 -.500000 E+00 

52 -.480000 E+00 

53 -.460000 E+00 

54 -440000 E+00 

55 -.420000 E+00 

56 -.400000 E+00 

57 -.380000 E+00 

58 -.360000 E+00 

59 -.340000 E+00 

60 -.320000 E+00 

61 -.300000 E+00 

62 -.280000 E+00 

63 -.260000 E+00 

64 -.240000 E+00 

65 -.220000 E+00 

66 -.200000 E+00 

67 -.180000 E+00 

68 -.160000 E+00 

69 -.140000 E+00 

70 -.120000 E+00 

71 -.100000 E+00 

72 -.800000 E-01 

73 -.600000 E-01 

74 -.400000 E-01 

75 -.200000 E-01 

76 0. 

77 .200000 E-01 

78 .400000 E-01 

79 .600000 E-01 

80 .800000 E-01 

81 .100000 E+00 

82 .120000 E+00 

83 .140000 E+00 

84 .160000 E+00 

85 .180000 E+00 

86 .200000 E+00 

87 .220000 E+00 

88 .240000 E+00 

89 .260000 E+00 

90 .280000 E+00 

91 .300000 E+00 

92 .320000 E+00 

93 .340000 E+00 

94 .360000 E+00 

95 .380000 E+00 

96 .400001 E+00 

97 .420000 E+00 

98 .440000 E+00 

99 .460000 E+00 

100 .480000 E+00 


i 

1 

101 

.500000 E+00 

102 

.520000 E+00 

103 

.540000 E+00 

104 

.560000 E+00 

105 

.580000 E+00 

106 

.600000 E+00 

107 

.620000 E+00 

108 

.640000 E+00 

109 

.660000 E+00 

110 

.680000 E+00 

111 

.700000 E+00 

112 

.720000 E+00 

1 13 

.740000 E+00 

1 14 

.760000 E+00 

115 

.780000 E+00 

116 

.800000 E+00 

117 

.820000 E+00 

118 

.840000 E+00 

119 

.860000 E+00 

120 

.880000 E+00 

121 

.900000 E+00 

122 

.920000 E+00 

123 

.940000 E+00 

124 

.960000 E+00 

125 

.980000 E+00 

126 

.100000 E+01 

127 

.102020 E+01 

128 

.104167 E+01 

129 

.106601 E+01 

130 

.109522 E+01 

131 

.113186 E+01 

132 

.117913 E+01 

133 

.124112 E+01 

134 

.132297 E+01 

135 

.143122 E+01 

136 

.157411 E+01 

137 

.176209 E+01 

138 

.200837 E+01 

139 

.232969 E+01 

140 

.274720 E+01 

141 

.328770 E+01 

142 

.398506 E+01 

143 

.488218 E+01 

144 

.603333 E+01 

145 

.750718 E+01 

146 

.939063 E+01 

147 

.1 17937 E+02 

148 

.148554 E+02 

149 

.187520 E+02 

150 

.237062 E+02 

151 

.300000 E+02 
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b) n ~ Direction Mesh Distribution 


V = 

y -S(x,t) 


j 

y 


1 

-.283650 

E-14 

2 

.188308 

E-01 

3 

.412027 

E-01 

4 

.677793 

E-01 

5 

.993481 

E-01 

6 

.136843 

E+00 

7 

.181370 

E+00 

8 

.234240 

E+00 

9 

.297005 

E+00 

10 

.371501 

E+00 

1 1 

.459898 

E+00 

12 

.564756 

E+00 

13 

.689097 

E+00 

14 

.836476 

E+00 

15 

.101107 

E+01 

16 

.121779 

E+01 

17 

.146236 

E+01 

18 

.175148 

E+01 

19 

.209290 

E+01 

20 

.249562 

E+01 

21 

.296999 

E+01 

22 

.352782 

E+01 

23 

.418255 

E+01 

24 

.494926 

E+01 

25 

.584473 

E+01 

26 

.688734 

E+01 

27 

.809692 

E+01 

28 

.949432 

E+01 

29 

.111010 

E+02 

30 

.129379 

E+02 

31 

.150249 

E+02 

32 

.173791 

E+02 

33 

.200129 

E+02 

34 

.229332 

E+02 

35 

.261386 

E+02 

36 

.296182 

E+02 

37 

.333506 

E+02 

38 

.373030 

E+02 

39 

.414318 

E+02 

40 

.456841 

E+02 

41 

.500000 

E+02 
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APPENDIX C 


EFFECTS OF VARIOUS APPROXIMATIONS TO FULL 
POTENTIAL THEORY ON PRESSURE COEFFICIENTS 


We have made additional computations of the flow field about a 
thickening and subsequently thinning airfoil. These computations use 
two modifications to full potential theory which are intermediate in 
going to the complete approximations made by low frequency, small dis- 
turbance, transonic theory. The results of the two modifications are 
shown in Figures 8 and 9 for two different times of the flow. 


The modification 1 results were obtained by replacing the actual 
airfoil shape boundary y = S(x,t) by the line segment y = 0, for 
0 X ^ 1. On that line segment, the airfoil boundary condition, Eq. (42) 
was replaced by the equation, 




y 


•f S ({) 
x^x 


(Cl) 


Modification 2 made the following assumptions. It also replaced the 
airfoil shape by the line segment y = 0. The airfoil boundary condition 
here became 


^ + S . (C2) 

y t X 

Also the formula for Cp , Eq. (16), was replaced by the low frequency, 
small disturbance, transonic formula 


Cp = -2 ((f)^ - 1) ; 

here (1)^ is the full potential and not the perturbation potential. Both 
modifications used the full potential equation, Eq. (18), as the governing 
equation. 

In Figure 8, we see that both modifications give shock locations in 
good agreement with full potential theory. However Modification 2, with 
its small disturbance boundary conditions, yields the least expansion 
upstream of this shock and the weakest shock. In those two aspects, the 
intermediate theory. Modification 2, has less agreement with full potential 
theory than complete small disturbance theory has. Modification 1 only 
lacks the physical displacement of the fluid caused by the thickness of 
the airfoil. Modification 1 has good agreement with full potential theory, 
although on the first 60% of the airfoil Modification 1 gives slightly 
more negative values of Cp. 
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■ ■ p ^ Full Potential Equation 

— * “ “ Modification 1 of 

Full Potential Equation 

““ Modification 2 of 

Full Potential Equation 

-**•• ^** Low Frequency, Small 
Disturbance, Transonic 
Equation, Ref. 6 


M^ * 0.85 
r* 18.25 


Figure 8.— Plot of Pressure Coefficients for Various Approximations to 
Full Potential Theory. Strong Shock Case. 


In Figure 9> the modifications are plotted at a later time when the 
shocks are now weak. Again, Modification 1 is in close agreement with 
full potential theory; now the airfoil has almost no thickness. Modifi- 
cation 2 again has less agreement with full potential theory in the 
shock region than the complete small disturbance approximation has. 
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- 0.1 - 


Full Potential Equation 

^ ^ ^ Modification 1 of 

Full Potential Equation 

Modification 2 of 

2 Full Potential Equation 

_■ _ Low Frequency, Small 
S.D. Disturbance, Transonic 
Equation, Ref. 6 


0.1 



26.875 


0.2 


- 1.0 


-0.5 


0 ^ 0.6 
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Figure 9.— Plot of Pressure Coefficients for Various Approximations to 
Full Potential Theory. Weak Shock Case. 
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